Here we are generating random steps for fitting step selection models, and plotting the outputs of the step generation. We are using a Gamma distribution for the step lengths and a von Mises distribution for the turning angles.
library(tidyverse)
packages <- c("amt", "lubridate", "terra", "beepr", "tictoc")
walk(packages, require, character.only = T)
## # A tibble: 6 × 4
## id time lon lat
## <chr> <dttm> <dbl> <dbl>
## 1 2005 2018-07-25 10:04:02 134. -12.3
## 2 2005 2018-07-25 11:04:23 134. -12.3
## 3 2005 2018-07-25 12:04:39 134. -12.3
## 4 2005 2018-07-25 13:04:17 134. -12.3
## 5 2005 2018-07-25 14:04:39 134. -12.3
## 6 2005 2018-07-25 15:04:27 134. -12.3
## [1] "Australia/Queensland"
Use the amt package to create a trajectory object from
the cleaned data.
buffalo_all <- buffalo_clean %>% mk_track(id = id,
lon,
lat,
time,
all_cols = T,
crs = 4326) %>%
transform_coords(crs_to = 3112, crs_from = 4326) # Transformation to GDA94 /
# Geoscience Australia Lambert (https://epsg.io/3112)
# plot the data spatially coloured by time
buffalo_all %>%
ggplot(aes(x = x_, y = y_, colour = t_)) +
geom_point(alpha = 0.25, size = 1) + # colour = "black",
coord_fixed() +
scale_colour_viridis_c() +
theme_classic()
To determine required extent of raster to crop environmental variables to (with a buffer to accommodate the random steps)
min(buffalo_all$x_)
## [1] 5337.494
max(buffalo_all$x_)
## [1] 52901.66
min(buffalo_all$y_)
## [1] -1457924
max(buffalo_all$y_)
## [1] -1411517
template_raster_cropped <- terra::rast(xmin = 0,
xmax = 60000,
ymin = -1463000,
ymax = -1406000,
resolution = 25,
crs = "epsg:3112")
template_wgs84 <- project(template_raster_cropped, "epsg:4326")
terra::ext(template_wgs84)
## SpatExtent : 134.000008561276, 134.541785217172, -12.5424223864104, -12.0361830606666 (xmin, xmax, ymin, ymax)
ndvi_projected <-
rast("mapping/cropped rasters/ndvi_GEE_projected_watermask20230207.tif")
slope <- rast("mapping/cropped rasters/slope_raster.tif")
veg_herby <- rast("mapping/cropped rasters/veg_herby.tif")
canopy_cover <- rast("mapping/cropped rasters/canopy_cover.tif")
# change the names (these will become the column names when extracting
# covariate values at the used and random steps)
names(ndvi_projected) <- rep("ndvi", terra::nlyr(ndvi_projected))
names(slope) <- "slope"
names(veg_herby) <- "veg_herby"
names(canopy_cover) <- "canopy_cover"
# plot the rasters
plot(ndvi_projected)
plot(slope)
plot(veg_herby)
plot(canopy_cover)
# nest the data by individual
buffalo_all_nested <- buffalo_all %>% arrange(id) %>% nest(data = -"id")
buffalo_all_nested_steps_by_burst <- buffalo_all_nested %>%
mutate(steps = map(data, function(x)
x %>% track_resample(rate = hours(1), tolerance = minutes(10)) %>%
# to filter out bursts with less than 3 locations
# amt::filter_min_n_burst(min_n = 3) %>%
steps_by_burst()))
# unnest the data after creating 'steps' objects
buffalo_all_steps_by_burst <- buffalo_all_nested_steps_by_burst %>%
amt::select(id, steps) %>%
amt::unnest(cols = steps)
buffalo_all_steps_by_burst <- buffalo_all_steps_by_burst %>%
mutate(t2_rounded = round_date(t2_, "hour"), # round the time to the nearest hour
hour_t2 = ifelse(hour(t2_rounded) == 0, 24, hour(t2_rounded))) # change the 0 hour to 24
head(buffalo_all_steps_by_burst, 10)
## # A tibble: 10 × 14
## id burst_ x1_ x2_ y1_ y2_ sl_ direction_p ta_ t1_ t2_ dt_
## * <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dttm> <dttm> <drt>
## 1 2005 1 41940. 41968. -1.44e6 -1.44e6 206. 1.43 NA 2018-07-25 10:04:02 2018-07-25 11:04:23 60.3…
## 2 2005 1 41968. 41921. -1.44e6 -1.44e6 50.7 2.80 1.37 2018-07-25 11:04:23 2018-07-25 12:04:39 60.2…
## 3 2005 1 41921. 41778. -1.44e6 -1.44e6 152. 2.78 -0.0214 2018-07-25 12:04:39 2018-07-25 13:04:17 59.6…
## 4 2005 1 41778. 41840. -1.44e6 -1.44e6 70.7 -0.507 2.99 2018-07-25 13:04:17 2018-07-25 14:04:39 60.3…
## 5 2005 1 41840. 41655. -1.44e6 -1.44e6 188. 2.98 -2.80 2018-07-25 14:04:39 2018-07-25 15:04:27 59.8…
## 6 2005 1 41655. 41618. -1.44e6 -1.44e6 37.1 -3.02 0.285 2018-07-25 15:04:27 2018-07-25 16:04:24 59.9…
## 7 2005 1 41618. 41687. -1.44e6 -1.44e6 522. -1.44 1.58 2018-07-25 16:04:24 2018-07-25 17:04:28 60.0…
## 8 2005 1 41687. 41683. -1.44e6 -1.44e6 8.13 2.10 -2.75 2018-07-25 17:04:28 2018-07-25 18:04:31 60.0…
## 9 2005 1 41683. 41674. -1.44e6 -1.44e6 9.75 -3.13 1.05 2018-07-25 18:04:31 2018-07-25 19:04:33 60.0…
## 10 2005 1 41674. 41424. -1.44e6 -1.44e6 308. 2.52 -0.635 2018-07-25 19:04:33 2018-07-25 20:04:02 59.4…
## # ℹ 2 more variables: t2_rounded <dttm>, hour_t2 <dbl>
Fitting exponential and von Mises distributions to the steps of ALL individuals (only one Gamma and one von Mises distribution for the whole population). This should be done when fitting a hierarchical model to update the ‘population’ parameters, but also makes it straightforward to update after model fitting to each individual separately.
# fitting step length and turning angle distributions to all locations
gamma_dist <- fit_distr(buffalo_all_steps_by_burst$sl_, "gamma")
vonmises_dist <- fit_distr(buffalo_all_steps_by_burst$ta_, "vonmises")
# checking parameters - which can then be saved to update movement parameters
# after fitting the step selection model
gamma_dist$params$shape
## [1] 0.438167
gamma_dist$params$scale
## [1] 534.3507
vonmises_dist$params$kappa
## [1] 0.1848126
vonmises_dist$params$mu
## Circular Data:
## Type = angles
## Units = radians
## Template = none
## Modulo = asis
## Zero = 0
## Rotation = counter
## [1] 0
For some reason, the random_steps function does not work
when using the bursted_steps_xyt class. I’m not sure why
(the error is
Error in bursts[[i]] : subscript out of bounds), but it
works when that class label is removed, and appears to sample random
steps correctly. For taxa that have irregular fixes and many bursts,
this may be worth exploring in more detail.
class(buffalo_all_steps_by_burst)
## [1] "bursted_steps_xyt" "steps_xyt" "steps_xy" "tbl_df" "tbl"
## [6] "data.frame"
class(buffalo_all_steps_by_burst) <- class(buffalo_all_steps_by_burst)[-1]
class(buffalo_all_steps_by_burst)
## [1] "steps_xyt" "steps_xy" "tbl_df" "tbl" "data.frame"
tic()
buffalo_parametric_popn_GvM <- buffalo_all_steps_by_burst %>%
random_steps(n_control = 10,
sl_distr = gamma_dist,
ta_distr = vonmises_dist) %>%
mutate(y = as.numeric(case_))
toc()
## 57.25 sec elapsed
Spatially plot the used and random steps, with the used steps in red.
buffalo_parametric_popn_GvM %>% ggplot() +
geom_point(data = . %>% filter(y == 0), aes(x = x2_, y = y2_),
colour = "black", size = 0.1, alpha = 0.1) +
geom_point(data = . %>% filter(y == 1), aes(x = x2_, y = y2_),
colour = "red", size = 0.1, alpha = 0.1) +
coord_equal() +
theme_bw()
Gamma distribution of the step lengths
buffalo_parametric_popn_GvM %>% ggplot() +
geom_density(data = . %>% filter(y == 1), aes(x = sl_, fill = factor(id)),
alpha = 0.25) +
geom_density(data = . %>% filter(y == 0), aes(x = sl_), colour = "red") +
scale_y_continuous("Density") +
scale_x_continuous("Step length (m)") +
scale_fill_viridis_d(direction = -1) +
theme_classic() +
theme(legend.position = "none")
ggsave(paste0("outputs/plots/manuscript_figs/step_length_density_", Sys.Date(), ".png"),
width=150, height=90, units="mm", dpi = 600)
When log-transforming the step length, it is clear that there is a bimodal pattern to the step lengths, which is not captured by the Gamma distribution of the random steps.
buffalo_parametric_popn_GvM %>% ggplot() +
geom_density(data = . %>% filter(y == 1), aes(x = log(sl_), fill = factor(id)),
alpha = 0.25) +
geom_density(data = . %>% filter(y == 0), aes(x = log(sl_)), colour = "red") +
scale_y_continuous("Density") +
scale_x_continuous("Log of step length (m)") +
scale_fill_viridis_d(direction = -1) +
theme_classic() +
theme(legend.position = "none")
ggsave(paste0("outputs/plots/manuscript_figs/step_length_density_log_", Sys.Date(), ".png"),
width=150, height=90, units="mm", dpi = 600)
Plot the von Mises distribution of the turning angles
buffalo_parametric_popn_GvM %>% ggplot() +
geom_density(data = . %>% filter(y == 1), aes(x = ta_, fill = factor(id)),
alpha = 0.25) +
geom_density(data = . %>% filter(y == 0), aes(x = ta_), colour = "red",
alpha = 0.01) +
scale_y_continuous("Density") +
scale_x_continuous("Turning angle (rad)") +
scale_fill_viridis_d(direction = -1) +
theme_classic() +
theme(legend.position = "none")
## Warning: Removed 8053 rows containing non-finite values (`stat_density()`).
ggsave(paste0("outputs/plots/manuscript_figs/turning angle_density_", Sys.Date(), ".png"),
width=150, height=90, units="mm", dpi = 600)
## Warning: Removed 8053 rows containing non-finite values (`stat_density()`).
buffalo_parametric_popn_covs <- buffalo_parametric_popn_GvM %>%
extract_covariates_var_time(ndvi_projected,
where = "end",
when = "any",
max_time = days(15),
name_covar = "ndvi_temporal") %>%
extract_covariates(veg_herby,
where = "end") %>%
extract_covariates(canopy_cover,
where = "end") %>%
extract_covariates(slope,
where = "end") %>%
mutate(y = as.numeric(case_),
cos_ta_ = cos(ta_),
log_sl_ = log(sl_))
# write the output to a csv file
write_csv(buffalo_parametric_popn_covs,
paste0("outputs/buffalo_parametric_popn_covs_GvM_10rs_", Sys.Date(), ".csv"))
beep(sound = 2)